###replication code for: political competition and repression####
#world politics
#pearce edwards
#april 25, 2022

#1. preparation####

rm(list = ls())
cat("\014")
set.seed(30030)

library(plm)
library(glm2)
library(reshape)
library(data.table)
library(ggthemes)
library(survival)
library(lmtest)
library(rgdal)
library(rgeos)
library(maptools)
library(ggplot2)
library(dplyr)
library(plyr)
library(stringr)
library(stringi)
library(gridExtra)
library(estimatr)
library(magrittr)
library(MASS)
library(pscl)
library(sensemakr)
library(xtable)
library(texreg)
library(mediation)

#1.1 load and clean data####

chile.commune <- read.csv("chile.commune.csv")
chile.commune.late <- read.csv("chile.commune.0913-1231.csv")
chile.commune.0973 <- read.csv("chile.commune.0911-0930.csv")
chile.commune.1073 <- read.csv("chile.commune.1001-1031.csv")
chile.commune.1173 <- read.csv("chile.commune.1101-1130.csv")
chile.commune.1273 <- read.csv("chile.commune.1201-1231.csv")
chile <- read.csv("chile.victimlevel.csv")

#1.2 descriptive graphs####

desc1 <- ggplot(data=subset(chile,!is.na(chile$profession_new)),aes(as.character(profession_new))) + 
  geom_bar(colour="black",fill="gray60") +
  ylab("Victim Count") + 
  xlab("Profession") + 
  ylim(0,800) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

desc2 <- ggplot(data=subset(chile,!is.na(chile$militancia)),aes(as.character(militancia))) + 
  geom_bar(colour="black", fill="gray60") +
  ylab("Victim Count") + 
  xlab("Dissident Classification") +  
  ylim(0,800) + 
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

desc.plot <- grid.arrange(desc1, desc2, ncol=2)

#figure 1
ggsave("descriptives.jpg", plot = desc.plot,
       scale = 3, width = 900, height = 400, units = "px",
       dpi = 300, limitsize = TRUE)

#2 maps####

#shape data from informacion teritorial - biblioteca del congreso nacional de chile
chile.shape <- readOGR(dsn=path.expand("division_comunal.shp"), layer =  "division_comunal")
chile.map.data <- fortify(chile.shape, region="NOM_COM")

chile.shape2 <- chile.shape@data

chile.map.data$id <- stri_trans_general(chile.map.data$id,"Latin-ASCII")
chile.map.data$id <- tolower(chile.map.data$id)
chile.map.data$id <- as.factor(chile.map.data$id)

chile.map.data$group <- stri_trans_general(chile.map.data$group,"Latin-ASCII")
chile.map.data$group <- tolower(chile.map.data$group)
chile.map.data$group <- as.factor(chile.map.data$group)

colnames(chile.map.data)[colnames(chile.map.data) == "id"] <- "comuna"

chile.shape2$NOM_COM <- as.character(chile.shape2$NOM_COM)
chile.shape2$NOM_COM <- stri_trans_general(chile.shape2$NOM_COM,"Latin-ASCII")
chile.shape2$NOM_COM <- tolower(chile.shape2$NOM_COM)
chile.shape2$NOM_COM <- as.factor(chile.shape2$NOM_COM)

colnames(chile.shape2)[colnames(chile.shape2) == "NOM_COM"] <- "comuna" 
chile.shape2$comuna <- as.character(chile.shape2$comuna)

chile.shape2 <- left_join(chile.shape2, chile.commune, by = "comuna")

map.data <- merge(chile.map.data, chile.shape2, by = "comuna", all.x = TRUE)
map.data <- subset(map.data, map.data$comuna != "isla de pascua")

map <- ggplot() +
  geom_polygon(data=map.data, aes(y=lat, x=long, group=group, fill=as.factor(ifelse(totalvictims == 0 | is.na(totalvictims),0,1))), color="gray90")

map2 <- ggplot() +
  geom_polygon(data=map.data, aes(y=lat, x=long, group=group, fill=as.factor(ifelse(elec1973!=2 | is.na(elec1973),0,1))), color="gray90")

map <- map + coord_equal()

map2 <- map2 + coord_equal()

map <- map + theme(axis.text.x = element_blank(),
                   axis.text.y = element_blank(),
                   axis.ticks = element_blank(),
                   axis.title.x = element_blank(),
                   axis.title.y = element_blank(),
                   axis.line.x = element_blank(),
                   panel.background = element_blank(),
                   panel.border = element_blank(),
                   panel.grid.major = element_blank(),
                   panel.grid.minor = element_blank(),
                   plot.background = element_blank())

map2 <- map2 + theme(axis.text.x = element_blank(),
                     axis.text.y = element_blank(),
                     axis.ticks = element_blank(),
                     axis.title.x = element_blank(),
                     axis.title.y = element_blank(),
                     axis.line.x = element_blank(),
                     panel.background = element_blank(),
                     panel.border = element_blank(),
                     panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(),
                     plot.background = element_blank())

map <- map + scale_fill_manual(values=c("gray95","gray35"), labels=c("No Killings","Political Killings")) + theme(legend.title=element_blank()) + 
  ggtitle("Political Killings") + theme(plot.title = element_text(size = 16, face = "bold"))

map2 <- map2 + scale_fill_manual(values=c("gray95","gray35"), labels=c("Not Close Election","Close Election")) + theme(legend.title=element_blank()) + 
  ggtitle("Close Elections, 1973") + theme(plot.title = element_text(size = 16, face = "bold"))

#figure 2
ggsave("map-killings.jpg", plot = map,
       scale = 1, width = 4, height = 8, units = "in",
       dpi = 300, limitsize = TRUE)

ggsave("map-closeelection.jpg", plot = map2,
       scale = 1, width = 4, height = 8, units = "in",
       dpi = 300, limitsize = TRUE)

#3. estimation####

#3.1 ols  ####

#total victims outcome
c_model1a <- lm_robust(log(totalvictims+1) ~ close_election_3pt, data=chile.commune,se_type="HC2")

c_model1b <- lm_robust(log(totalvictims+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                           as.factor(province) + armybase, data=chile.commune, se_type="HC2")

#non-militant outcome
c_model1c <- lm_robust(log(nonmilitant+1) ~ close_election_3pt, data=chile.commune,se_type="HC2")

c_model1d <- lm_robust(log(nonmilitant+1) ~ close_election_3pt  + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")

#substantive significance
as.numeric(c_model1b$coefficients[2])/sd(log(chile.commune$totalvictims+1),na.rm=T)
as.numeric(c_model1d$coefficients[2])/sd(log(chile.commune$totalvictims+1),na.rm=T)

exp(c_model1b$coefficients[2])
exp(c_model1d$coefficients[2])

#plot
mass.plot2 <- matrix(NA,4,3)

mass.plot2[,1] <- c(as.numeric(c_model1a$coefficients[2]) - 1.96*as.numeric(c_model1a$std.error[2]),
                    as.numeric(c_model1b$coefficients[2]) - 1.96*as.numeric(c_model1b$std.error[2]),
                    as.numeric(c_model1c$coefficients[2]) - 1.96*as.numeric(c_model1c$std.error[2]),
                    as.numeric(c_model1d$coefficients[2]) - 1.96*as.numeric(c_model1d$std.error[2]))
mass.plot2[,2] <- c(as.numeric(c_model1a$coefficients[2]),
                    as.numeric(c_model1b$coefficients[2]),
                    as.numeric(c_model1c$coefficients[2]),
                    as.numeric(c_model1d$coefficients[2]))
mass.plot2[,3] <- c(as.numeric(c_model1a$coefficients[2]) + 1.96*as.numeric(c_model1a$std.error[2]),
                    as.numeric(c_model1b$coefficients[2]) + 1.96*as.numeric(c_model1b$std.error[2]),
                    as.numeric(c_model1c$coefficients[2]) + 1.96*as.numeric(c_model1c$std.error[2]),
                    as.numeric(c_model1d$coefficients[2]) + 1.96*as.numeric(c_model1d$std.error[2]))


mass.plot2 <- as.data.frame(mass.plot2)
mass.plot2$names <- c("Total Killings","Total Killings, Covariates","Non-Dissidents","Non-Dissidents, Covariates")
mass.plot2$names <- as.factor(mass.plot2$names)
mass.plot2$names <- factor(mass.plot2$names,levels(mass.plot2$names)[c(2,1,4,3)])
colnames(mass.plot2) <- c("lower","ptest","upper","names")

mass.cfplot2<- ggplot(data = mass.plot2, aes(x = ptest, y = names, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
  geom_vline(xintercept =0, size=1, linetype = "dashed") +
  xlim(-0.25,1) + 
  xlab("Estimated Change in Log(Political Killings)") +
  ylab("Dependent Variable") +
  theme_few() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("black","gray35","black","gray35"))

mass.cfplot2

#figure 3
ggsave("plot-closecomp-1.jpg", plot = mass.cfplot2,
       scale = 3, width = 550, height = 504, units = "px",
       dpi = 300, limitsize = TRUE)

#3.2 robustness: count models ####

#negative binomial ####

c_model3a <- glm.nb(totalvictims ~ close_election_3pt, data=chile.commune)

c_model3b <- glm.nb(totalvictims ~ close_election_3pt + allende_1970_pct + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                      as.factor(province) + armybase, data=chile.commune)


c_model3c <- glm.nb(nonmilitant ~ close_election_3pt, data=chile.commune)

c_model3d <- glm.nb(nonmilitant ~ close_election_3pt + allende_1970_pct + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                      as.factor(province) + armybase, data=chile.commune)

mass.plot <- matrix(NA,4,3)

mass.plot[,1] <- c(summary(c_model3a)$coefficients[2,1] - 1.96*summary(c_model3a)$coefficients[2,2],
                   summary(c_model3b)$coefficients[2,1] - 1.96*summary(c_model3b)$coefficients[2,2], 
                   summary(c_model3c)$coefficients[2,1] - 1.96*summary(c_model3c)$coefficients[2,2], 
                   summary(c_model3d)$coefficients[2,1] - 1.96*summary(c_model3d)$coefficients[2,2])
mass.plot[,2] <- c(summary(c_model3a)$coefficients[2,1],
                   summary(c_model3b)$coefficients[2,1],
                   summary(c_model3c)$coefficients[2,1],
                   summary(c_model3d)$coefficients[2,1])
mass.plot[,3] <- c(summary(c_model3a)$coefficients[2,1] + 1.96*summary(c_model3a)$coefficients[2,2],
                   summary(c_model3b)$coefficients[2,1] + 1.96*summary(c_model3b)$coefficients[2,2], 
                   summary(c_model3c)$coefficients[2,1] + 1.96*summary(c_model3c)$coefficients[2,2], 
                   summary(c_model3d)$coefficients[2,1] + 1.96*summary(c_model3d)$coefficients[2,2])

mass.plot <- as.data.frame(mass.plot)
mass.plot$names <- c("Total Killings","Total Killings, Covariates","Non-Dissidents","Non-Dissidents, Covariates")
mass.plot$names <- as.factor(mass.plot$names)
mass.plot$names <- factor(mass.plot$names,levels(mass.plot$names)[c(2,1,4,3)])
colnames(mass.plot) <- c("lower","ptest","upper","names")

mass.cfplot<- ggplot(data = mass.plot, aes(x = ptest, y = names, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(x = lower, xend = upper, y = names, yend = names), size=1.25) +
  geom_vline(xintercept =0, size=1, linetype = "dashed") +
  xlim(-1,3) + 
  xlab("Estimated Change in Political Killings") +
  ylab("Dependent Variable") +
  theme_few() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("black","gray35","black","gray35"))

mass.cfplot

#zero-inflated negative binomial ####

c_model4a <- zeroinfl(totalvictims ~ close_election_3pt + log(pop_1970_est) | log(pop_1970_est), data=chile.commune, dist = "negbin")

c_model4b <- zeroinfl(totalvictims ~ close_election_3pt + allende_1970_pct + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                        as.factor(province) + armybase | log(pop_1970_est), data=chile.commune, dist = "negbin")

c_model4c <- zeroinfl(nonmilitant ~ close_election_3pt + log(pop_1970_est) | log(pop_1970_est), data=chile.commune, dist = "negbin")

c_model4d <- zeroinfl(nonmilitant ~ close_election_3pt + allende_1970_votes + polarized + infant_mortality_pct  + log(pop_1970_est) + votershare +
                        as.factor(province) + armybase | log(pop_1970_est), data=chile.commune, dist = "negbin")

#3.3 robustness: dv killings per 100k ####

c_model5a <- lm_robust(log(totalvictims_per100k+1) ~ close_election_3pt, data=chile.commune,se_type="HC2")

c_model5b <- lm_robust(log(totalvictims_per100k+1) ~ close_election_3pt + allende_1970_pct + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")

c_model5c <- lm_robust(log(nonmilitant_per100k+1) ~ close_election_3pt, data=chile.commune,se_type="HC2")

c_model5d <- lm_robust(log(nonmilitant_per100k+1) ~ close_election_3pt + allende_1970_pct + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")

#3.4 robustness: timing ####

#sept 13 or later ####

#with up 1973
c_model13a <- lm_robust(log(totalvictims+1) ~ close_election_3pt, data=chile.commune.late,se_type="HC2")

c_model13b <- lm_robust(log(totalvictims+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune.late, se_type="HC2")

c_model13c <- lm_robust(log(nonmilitant+1) ~ close_election_3pt, data=chile.commune.late,se_type="HC2")

c_model13d <- lm_robust(log(nonmilitant+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune.late, se_type="HC2")


#by-month effects ####

#total victims outcome

#september
c_model13e <- lm_robust(log(totalvictims+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune.0973, se_type="HC2")

#october
c_model13f <- lm_robust(log(totalvictims+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune.1073, se_type="HC2")

#november
c_model13g <- lm_robust(log(totalvictims+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune.1173, se_type="HC2")

#december
c_model13h <- lm_robust(log(totalvictims+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune.1273, se_type="HC2")

#non-militant outcome

#september
c_model13i <- lm_robust(log(nonmilitant+1) ~ close_election_3pt  + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune.0973, se_type="HC2")

#october
c_model13j <- lm_robust(log(nonmilitant+1) ~ close_election_3pt  + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune.1073, se_type="HC2")

#november
c_model13k <- lm_robust(log(nonmilitant+1) ~ close_election_3pt  + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune.1173, se_type="HC2")

#december
c_model13l <- lm_robust(log(nonmilitant+1) ~ close_election_3pt  + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune.1273, se_type="HC2")

#plot
mass.plot3 <- matrix(NA,4,3)

mass.plot3[,1] <- c(summary(c_model13e)$coefficients[2,1] - 1.96*summary(c_model13e)$coefficients[2,2],
                    summary(c_model13f)$coefficients[2,1] - 1.96*summary(c_model13f)$coefficients[2,2], 
                    summary(c_model13g)$coefficients[2,1] - 1.96*summary(c_model13g)$coefficients[2,2], 
                    summary(c_model13h)$coefficients[2,1] - 1.96*summary(c_model13h)$coefficients[2,2])
mass.plot3[,2] <- c(summary(c_model13e)$coefficients[2,1],
                    summary(c_model13f)$coefficients[2,1], 
                    summary(c_model13g)$coefficients[2,1], 
                    summary(c_model13h)$coefficients[2,1])
mass.plot3[,3] <- c(summary(c_model13e)$coefficients[2,1] + 1.96*summary(c_model13e)$coefficients[2,2],
                    summary(c_model13f)$coefficients[2,1] + 1.96*summary(c_model13f)$coefficients[2,2], 
                    summary(c_model13g)$coefficients[2,1] + 1.96*summary(c_model13g)$coefficients[2,2], 
                    summary(c_model13h)$coefficients[2,1] + 1.96*summary(c_model13h)$coefficients[2,2])

mass.plot3 <- as.data.frame(mass.plot3)
mass.plot3$names <- c("September","October","November","December")
mass.plot3$names <- as.factor(mass.plot3$names)
mass.plot3$names <- factor(mass.plot3$names,levels(mass.plot3$names)[c(4:1)])
colnames(mass.plot3) <- c("lower","ptest","upper","names")

mass.cfplot3 <- ggplot(data = mass.plot3, aes(y = ptest, x = names, color = names)) +
  geom_point(size=3.5) +
  geom_segment(aes(y = lower, yend = upper, x = names, xend = names), size=1.5) +
  geom_hline(yintercept =0, size=1, linetype = "dashed") +
  #xlim(-3,5.5) + 
  ylab("Estimated Change in Log(Political Killings)") +
  xlab("Month of 1973") +
  theme_few() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, vjust = 1.1, hjust=1)) +
  scale_color_manual(values=c("black","black","black","black","black","black","black","black","black","black",
                              "black", "black", "black", "black", "black", "black"))

mass.cfplot3

#figure 6
ggsave("plot-months-1.jpg", plot = mass.cfplot3,
       scale = 3, width = 550, height = 504, units = "px",
       dpi = 300, limitsize = TRUE)

mass.plot3b <- matrix(NA,4,3)

mass.plot3b[,1] <- c(summary(c_model13i)$coefficients[2,1] - 1.96*summary(c_model13i)$coefficients[2,2],
                    summary(c_model13j)$coefficients[2,1] - 1.96*summary(c_model13j)$coefficients[2,2], 
                    summary(c_model13k)$coefficients[2,1] - 1.96*summary(c_model13k)$coefficients[2,2], 
                    summary(c_model13l)$coefficients[2,1] - 1.96*summary(c_model13l)$coefficients[2,2])
mass.plot3b[,2] <- c(summary(c_model13i)$coefficients[2,1],
                    summary(c_model13j)$coefficients[2,1], 
                    summary(c_model13k)$coefficients[2,1], 
                    summary(c_model13l)$coefficients[2,1])
mass.plot3b[,3] <- c(summary(c_model13i)$coefficients[2,1] + 1.96*summary(c_model13i)$coefficients[2,2],
                    summary(c_model13j)$coefficients[2,1] + 1.96*summary(c_model13j)$coefficients[2,2], 
                    summary(c_model13k)$coefficients[2,1] + 1.96*summary(c_model13k)$coefficients[2,2], 
                    summary(c_model13l)$coefficients[2,1] + 1.96*summary(c_model13l)$coefficients[2,2])

mass.plot3b <- as.data.frame(mass.plot3b)
mass.plot3b$names <- c("September","October","November","December")
mass.plot3b$names <- as.factor(mass.plot3b$names)
mass.plot3b$names <- factor(mass.plot3b$names,levels(mass.plot3b$names)[c(4:1)])
colnames(mass.plot3b) <- c("lower","ptest","upper","names")

mass.cfplot3b <- ggplot(data = mass.plot3b, aes(y = ptest, x = names, color = names)) +
  geom_point(size=3.5) +
  geom_segment(aes(y = lower, yend = upper, x = names, xend = names), size=1.5) +
  geom_hline(yintercept =0, size=1, linetype = "dashed") +
  #xlim(-3,5.5) + 
  ylab("Estimated Change in Log(Non-Dissident Victims)") +
  xlab("Month of 1973") +
  theme_few() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, vjust = 1.1, hjust=1)) +
  scale_color_manual(values=c("black","black","black","black","black","black","black","black","black","black",
                              "black", "black", "black", "black", "black", "black"))

mass.cfplot3b

#figure 6
ggsave("plot-months-2.jpg", plot = mass.cfplot3b,
       scale = 3, width = 550, height = 504, units = "px",
       dpi = 300, limitsize = TRUE)

#3.5 robustness: placebo close election thresholds ####

#widening window ####

#total victims
c_model20 <- lm_robust(log(totalvictims+1) ~ close_election_17pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model21 <- lm_robust(log(totalvictims+1) ~ close_election_16pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model22 <- lm_robust(log(totalvictims+1) ~ close_election_15pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model23 <- lm_robust(log(totalvictims+1) ~ close_election_14pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model24 <- lm_robust(log(totalvictims+1) ~ close_election_13pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model25 <- lm_robust(log(totalvictims+1) ~ close_election_12pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model26 <- lm_robust(log(totalvictims+1) ~ close_election_11pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model27 <- lm_robust(log(totalvictims+1) ~ close_election_10pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model28 <- lm_robust(log(totalvictims+1) ~ close_election_9pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2") 
c_model29 <- lm_robust(log(totalvictims+1) ~ close_election_8pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model30 <- lm_robust(log(totalvictims+1) ~ close_election_7pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model31 <- lm_robust(log(totalvictims+1) ~ close_election_6pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model32 <- lm_robust(log(totalvictims+1) ~ close_election_5pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model33 <- lm_robust(log(totalvictims+1) ~ close_election_4pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model34 <- lm_robust(log(totalvictims+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model35 <- lm_robust(log(totalvictims+1) ~ close_election_2pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                  as.factor(province) + armybase, data=chile.commune, se_type="HC2")

mass.plot4 <- matrix(NA,16,3)

mass.plot4[,1] <- c(summary(c_model20)$coefficients[2,1] - 1.96*summary(c_model20)$coefficients[2,2],
                   summary(c_model21)$coefficients[2,1] - 1.96*summary(c_model21)$coefficients[2,2], 
                   summary(c_model22)$coefficients[2,1] - 1.96*summary(c_model22)$coefficients[2,2], 
                   summary(c_model23)$coefficients[2,1] - 1.96*summary(c_model23)$coefficients[2,2],
                   summary(c_model24)$coefficients[2,1] - 1.96*summary(c_model24)$coefficients[2,2],
                   summary(c_model25)$coefficients[2,1] - 1.96*summary(c_model25)$coefficients[2,2],
                   summary(c_model26)$coefficients[2,1] - 1.96*summary(c_model26)$coefficients[2,2],
                   summary(c_model27)$coefficients[2,1] - 1.96*summary(c_model27)$coefficients[2,2],
                   summary(c_model28)$coefficients[2,1] - 1.96*summary(c_model28)$coefficients[2,2],
                   summary(c_model29)$coefficients[2,1] - 1.96*summary(c_model29)$coefficients[2,2],
                   summary(c_model30)$coefficients[2,1] - 1.96*summary(c_model30)$coefficients[2,2],
                   summary(c_model31)$coefficients[2,1] - 1.96*summary(c_model31)$coefficients[2,2],
                   summary(c_model32)$coefficients[2,1] - 1.96*summary(c_model32)$coefficients[2,2],
                   summary(c_model33)$coefficients[2,1] - 1.96*summary(c_model33)$coefficients[2,2],
                   summary(c_model34)$coefficients[2,1] - 1.96*summary(c_model34)$coefficients[2,2],
                   summary(c_model35)$coefficients[2,1] - 1.96*summary(c_model35)$coefficients[2,2])
mass.plot4[,2] <- c(summary(c_model20)$coefficients[2,1],
                    summary(c_model21)$coefficients[2,1], 
                    summary(c_model22)$coefficients[2,1], 
                    summary(c_model23)$coefficients[2,1],
                    summary(c_model24)$coefficients[2,1],
                    summary(c_model25)$coefficients[2,1],
                    summary(c_model26)$coefficients[2,1],
                    summary(c_model27)$coefficients[2,1],
                    summary(c_model28)$coefficients[2,1],
                    summary(c_model29)$coefficients[2,1],
                    summary(c_model30)$coefficients[2,1],
                    summary(c_model31)$coefficients[2,1],
                    summary(c_model32)$coefficients[2,1],
                    summary(c_model33)$coefficients[2,1],
                    summary(c_model34)$coefficients[2,1],
                    summary(c_model35)$coefficients[2,1])
mass.plot4[,3] <- c(summary(c_model20)$coefficients[2,1] + 1.96*summary(c_model20)$coefficients[2,2],
                    summary(c_model21)$coefficients[2,1] + 1.96*summary(c_model21)$coefficients[2,2], 
                    summary(c_model22)$coefficients[2,1] + 1.96*summary(c_model22)$coefficients[2,2], 
                    summary(c_model23)$coefficients[2,1] + 1.96*summary(c_model23)$coefficients[2,2],
                    summary(c_model24)$coefficients[2,1] + 1.96*summary(c_model24)$coefficients[2,2],
                    summary(c_model25)$coefficients[2,1] + 1.96*summary(c_model25)$coefficients[2,2],
                    summary(c_model26)$coefficients[2,1] + 1.96*summary(c_model26)$coefficients[2,2],
                    summary(c_model27)$coefficients[2,1] + 1.96*summary(c_model27)$coefficients[2,2],
                    summary(c_model28)$coefficients[2,1] + 1.96*summary(c_model28)$coefficients[2,2],
                    summary(c_model29)$coefficients[2,1] + 1.96*summary(c_model29)$coefficients[2,2],
                    summary(c_model30)$coefficients[2,1] + 1.96*summary(c_model30)$coefficients[2,2],
                    summary(c_model31)$coefficients[2,1] + 1.96*summary(c_model31)$coefficients[2,2],
                    summary(c_model32)$coefficients[2,1] + 1.96*summary(c_model32)$coefficients[2,2],
                    summary(c_model33)$coefficients[2,1] + 1.96*summary(c_model33)$coefficients[2,2],
                    summary(c_model34)$coefficients[2,1] + 1.96*summary(c_model34)$coefficients[2,2],
                    summary(c_model35)$coefficients[2,1] + 1.96*summary(c_model35)$coefficients[2,2])

mass.plot4 <- as.data.frame(mass.plot4)
mass.plot4$names <- c("67-33","66-34","65-35","64-36","63-37","62-38","61-39","60-40","59-41", "58-42","57-43", "56-44", "55-45",
                      "54-46", "53-47", "52-48")
mass.plot4$names <- as.factor(mass.plot4$names)
#mass.plot4$names <- factor(mass.plot4$names,levels(mass.plot4$names)[c(9:16,1:8)])
colnames(mass.plot4) <- c("lower","ptest","upper","names")

mass.cfplot4 <- ggplot(data = mass.plot4, aes(y = ptest, x = names, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(y = lower, yend = upper, x = names, xend = names), size=0.5) +
  geom_hline(yintercept =0, size=1, linetype = "dashed") +
  #xlim(-3,5.5) + 
  ylab("Estimated Change in Political Killings") +
  xlab("UP Vote Share Range") +
  theme_few() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("black","black","black","black","black","black","black","black","black","black",
                              "black", "black", "black", "black", "black", "black"))

mass.cfplot4

#non-dissident victims
c_model20a <- lm_robust(log(nonmilitant+1) ~ close_election_17pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model21a <- lm_robust(log(nonmilitant+1) ~ close_election_16pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model22a <- lm_robust(log(nonmilitant+1) ~ close_election_15pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model23a <- lm_robust(log(nonmilitant+1) ~ close_election_14pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model24a <- lm_robust(log(nonmilitant+1) ~ close_election_13pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model25a <- lm_robust(log(nonmilitant+1) ~ close_election_12pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model26a <- lm_robust(log(nonmilitant+1) ~ close_election_11pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model27a <- lm_robust(log(nonmilitant+1) ~ close_election_10pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model28a <- lm_robust(log(nonmilitant+1) ~ close_election_9pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2") 
c_model29a <- lm_robust(log(nonmilitant+1) ~ close_election_8pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model30a <- lm_robust(log(nonmilitant+1) ~ close_election_7pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model31a <- lm_robust(log(nonmilitant+1) ~ close_election_6pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model32a <- lm_robust(log(nonmilitant+1) ~ close_election_5pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model33a <- lm_robust(log(nonmilitant+1) ~ close_election_4pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model34a <- lm_robust(log(nonmilitant+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model35a <- lm_robust(log(nonmilitant+1) ~ close_election_2pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")

mass.plot5 <- matrix(NA,16,3)

mass.plot5[,1] <- c(summary(c_model20a)$coefficients[2,1] - 1.96*summary(c_model20a)$coefficients[2,2],
                    summary(c_model21a)$coefficients[2,1] - 1.96*summary(c_model21a)$coefficients[2,2], 
                    summary(c_model22a)$coefficients[2,1] - 1.96*summary(c_model22a)$coefficients[2,2], 
                    summary(c_model23a)$coefficients[2,1] - 1.96*summary(c_model23a)$coefficients[2,2],
                    summary(c_model24a)$coefficients[2,1] - 1.96*summary(c_model24a)$coefficients[2,2],
                    summary(c_model25a)$coefficients[2,1] - 1.96*summary(c_model25a)$coefficients[2,2],
                    summary(c_model26a)$coefficients[2,1] - 1.96*summary(c_model26a)$coefficients[2,2],
                    summary(c_model27a)$coefficients[2,1] - 1.96*summary(c_model27a)$coefficients[2,2],
                    summary(c_model28a)$coefficients[2,1] - 1.96*summary(c_model28a)$coefficients[2,2],
                    summary(c_model29a)$coefficients[2,1] - 1.96*summary(c_model29a)$coefficients[2,2],
                    summary(c_model30a)$coefficients[2,1] - 1.96*summary(c_model30a)$coefficients[2,2],
                    summary(c_model31a)$coefficients[2,1] - 1.96*summary(c_model31a)$coefficients[2,2],
                    summary(c_model32a)$coefficients[2,1] - 1.96*summary(c_model32a)$coefficients[2,2],
                    summary(c_model33a)$coefficients[2,1] - 1.96*summary(c_model33a)$coefficients[2,2],
                    summary(c_model34a)$coefficients[2,1] - 1.96*summary(c_model34a)$coefficients[2,2],
                    summary(c_model35a)$coefficients[2,1] - 1.96*summary(c_model35a)$coefficients[2,2])
mass.plot5[,2] <- c(summary(c_model20a)$coefficients[2,1],
                    summary(c_model21a)$coefficients[2,1], 
                    summary(c_model22a)$coefficients[2,1], 
                    summary(c_model23a)$coefficients[2,1],
                    summary(c_model24a)$coefficients[2,1],
                    summary(c_model25a)$coefficients[2,1],
                    summary(c_model26a)$coefficients[2,1],
                    summary(c_model27a)$coefficients[2,1],
                    summary(c_model28a)$coefficients[2,1],
                    summary(c_model29a)$coefficients[2,1],
                    summary(c_model30a)$coefficients[2,1],
                    summary(c_model31a)$coefficients[2,1],
                    summary(c_model32a)$coefficients[2,1],
                    summary(c_model33a)$coefficients[2,1],
                    summary(c_model34a)$coefficients[2,1],
                    summary(c_model35a)$coefficients[2,1])
mass.plot5[,3] <- c(summary(c_model20a)$coefficients[2,1] + 1.96*summary(c_model20a)$coefficients[2,2],
                    summary(c_model21a)$coefficients[2,1] + 1.96*summary(c_model21a)$coefficients[2,2], 
                    summary(c_model22a)$coefficients[2,1] + 1.96*summary(c_model22a)$coefficients[2,2], 
                    summary(c_model23a)$coefficients[2,1] + 1.96*summary(c_model23a)$coefficients[2,2],
                    summary(c_model24a)$coefficients[2,1] + 1.96*summary(c_model24a)$coefficients[2,2],
                    summary(c_model25a)$coefficients[2,1] + 1.96*summary(c_model25a)$coefficients[2,2],
                    summary(c_model26a)$coefficients[2,1] + 1.96*summary(c_model26a)$coefficients[2,2],
                    summary(c_model27a)$coefficients[2,1] + 1.96*summary(c_model27a)$coefficients[2,2],
                    summary(c_model28a)$coefficients[2,1] + 1.96*summary(c_model28a)$coefficients[2,2],
                    summary(c_model29a)$coefficients[2,1] + 1.96*summary(c_model29a)$coefficients[2,2],
                    summary(c_model30a)$coefficients[2,1] + 1.96*summary(c_model30a)$coefficients[2,2],
                    summary(c_model31a)$coefficients[2,1] + 1.96*summary(c_model31a)$coefficients[2,2],
                    summary(c_model32a)$coefficients[2,1] + 1.96*summary(c_model32a)$coefficients[2,2],
                    summary(c_model33a)$coefficients[2,1] + 1.96*summary(c_model33a)$coefficients[2,2],
                    summary(c_model34a)$coefficients[2,1] + 1.96*summary(c_model34a)$coefficients[2,2],
                    summary(c_model35a)$coefficients[2,1] + 1.96*summary(c_model35a)$coefficients[2,2])

mass.plot5 <- as.data.frame(mass.plot5)
mass.plot5$names <- c("67-33","66-34","65-35","64-36","63-37","62-38","61-39","60-40","59-41", "58-42","57-43", "56-44", "55-45",
                      "54-46", "53-47", "52-48")
mass.plot5$names <- as.factor(mass.plot4$names)
#mass.plot5$names <- factor(mass.plot4$names,levels(mass.plot4$names)[c(9:16,1:8)])
colnames(mass.plot5) <- c("lower","ptest","upper","names")

mass.cfplot5 <- ggplot(data = mass.plot5, aes(y = ptest, x = names, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(y = lower, yend = upper, x = names, xend = names), size=0.5) +
  geom_hline(yintercept =0, size=1, linetype = "dashed") +
  #xlim(-3,5.5) + 
  ylab("Estimated Change in Non-Dissident Victims") +
  xlab("UP Vote Share Range") +
  theme_few() +
  theme(legend.position = "none") +
  scale_color_manual(values=c("black","black","black","black","black","black","black","black","black","black",
                              "black", "black", "black", "black", "black", "black"))

mass.cfplot5


#shifting windows ####

#total victims
c_model20c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= -5 & margin >= -11,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                            as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model21c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= -4 & margin >= -10,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                        as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model22c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= -3 & margin >= -9,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model23c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= -2 & margin >= -8,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model24c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= -1 & margin >= -7,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model25c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= 0 & margin >= -6,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model26c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= 1 & margin >= -5,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model27c <- lm_robust(log(totalvictims+1) ~  ifelse(margin <= 2 & margin >= -4,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model28c <- lm_robust(log(totalvictims+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model29c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= 4 & margin >= -2,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model30c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= 5 & margin >= -1,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model31c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= 6 & margin >= 0,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model32c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= 7 & margin >= 1,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model33c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= 8 & margin >= 2,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model34c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= 9 & margin >= 3,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model35c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= 10 & margin >= 4,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model36c <- lm_robust(log(totalvictims+1) ~ ifelse(margin <= 11 & margin >= 5,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")

mass.plot6 <- matrix(NA,17,3)

mass.plot6[,1] <- c(summary(c_model20c)$coefficients[2,1] - 1.96*summary(c_model20c)$coefficients[2,2],
                    summary(c_model21c)$coefficients[2,1] - 1.96*summary(c_model21c)$coefficients[2,2], 
                    summary(c_model22c)$coefficients[2,1] - 1.96*summary(c_model22c)$coefficients[2,2], 
                    summary(c_model23c)$coefficients[2,1] - 1.96*summary(c_model23c)$coefficients[2,2],
                    summary(c_model24c)$coefficients[2,1] - 1.96*summary(c_model24c)$coefficients[2,2],
                    summary(c_model25c)$coefficients[2,1] - 1.96*summary(c_model25c)$coefficients[2,2],
                    summary(c_model26c)$coefficients[2,1] - 1.96*summary(c_model26c)$coefficients[2,2],
                    summary(c_model27c)$coefficients[2,1] - 1.96*summary(c_model27c)$coefficients[2,2],
                    summary(c_model28c)$coefficients[2,1] - 1.96*summary(c_model28c)$coefficients[2,2],
                    summary(c_model29c)$coefficients[2,1] - 1.96*summary(c_model29c)$coefficients[2,2],
                    summary(c_model30c)$coefficients[2,1] - 1.96*summary(c_model30c)$coefficients[2,2],
                    summary(c_model31c)$coefficients[2,1] - 1.96*summary(c_model31c)$coefficients[2,2],
                    summary(c_model32c)$coefficients[2,1] - 1.96*summary(c_model32c)$coefficients[2,2],
                    summary(c_model33c)$coefficients[2,1] - 1.96*summary(c_model33c)$coefficients[2,2],
                    summary(c_model34c)$coefficients[2,1] - 1.96*summary(c_model34c)$coefficients[2,2],
                    summary(c_model35c)$coefficients[2,1] - 1.96*summary(c_model35c)$coefficients[2,2],
                    summary(c_model36c)$coefficients[2,1] - 1.96*summary(c_model36c)$coefficients[2,2])
mass.plot6[,2] <- c(summary(c_model20c)$coefficients[2,1],
                    summary(c_model21c)$coefficients[2,1], 
                    summary(c_model22c)$coefficients[2,1], 
                    summary(c_model23c)$coefficients[2,1],
                    summary(c_model24c)$coefficients[2,1],
                    summary(c_model25c)$coefficients[2,1],
                    summary(c_model26c)$coefficients[2,1],
                    summary(c_model27c)$coefficients[2,1],
                    summary(c_model28c)$coefficients[2,1],
                    summary(c_model29c)$coefficients[2,1],
                    summary(c_model30c)$coefficients[2,1],
                    summary(c_model31c)$coefficients[2,1],
                    summary(c_model32c)$coefficients[2,1],
                    summary(c_model33c)$coefficients[2,1],
                    summary(c_model34c)$coefficients[2,1],
                    summary(c_model35c)$coefficients[2,1],
                    summary(c_model36c)$coefficients[2,1])
mass.plot6[,3] <- c(summary(c_model20c)$coefficients[2,1] + 1.96*summary(c_model20c)$coefficients[2,2],
                    summary(c_model21c)$coefficients[2,1] + 1.96*summary(c_model21c)$coefficients[2,2], 
                    summary(c_model22c)$coefficients[2,1] + 1.96*summary(c_model22c)$coefficients[2,2], 
                    summary(c_model23c)$coefficients[2,1] + 1.96*summary(c_model23c)$coefficients[2,2],
                    summary(c_model24c)$coefficients[2,1] + 1.96*summary(c_model24c)$coefficients[2,2],
                    summary(c_model25c)$coefficients[2,1] + 1.96*summary(c_model25c)$coefficients[2,2],
                    summary(c_model26c)$coefficients[2,1] + 1.96*summary(c_model26c)$coefficients[2,2],
                    summary(c_model27c)$coefficients[2,1] + 1.96*summary(c_model27c)$coefficients[2,2],
                    summary(c_model28c)$coefficients[2,1] + 1.96*summary(c_model28c)$coefficients[2,2],
                    summary(c_model29c)$coefficients[2,1] + 1.96*summary(c_model29c)$coefficients[2,2],
                    summary(c_model30c)$coefficients[2,1] + 1.96*summary(c_model30c)$coefficients[2,2],
                    summary(c_model31c)$coefficients[2,1] + 1.96*summary(c_model31c)$coefficients[2,2],
                    summary(c_model32c)$coefficients[2,1] + 1.96*summary(c_model32c)$coefficients[2,2],
                    summary(c_model33c)$coefficients[2,1] + 1.96*summary(c_model33c)$coefficients[2,2],
                    summary(c_model34c)$coefficients[2,1] + 1.96*summary(c_model34c)$coefficients[2,2],
                    summary(c_model35c)$coefficients[2,1] + 1.96*summary(c_model35c)$coefficients[2,2],
                    summary(c_model36c)$coefficients[2,1] + 1.96*summary(c_model36c)$coefficients[2,2])

mass.plot6 <- as.data.frame(mass.plot6)
mass.plot6$names <- c("39-45","40-46","41-47","42-48","43-49","44-50","45-51","46-52","47-53","48-54",
                      "49-55","50-56","51-57","52-58","53-59","54-60","55-61")
mass.plot6$names <- as.factor(mass.plot6$names)
colnames(mass.plot6) <- c("lower","ptest","upper","names")

mass.cfplot6 <- ggplot(data = mass.plot6, aes(x = names, y = ptest, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(x = names, xend = names, y = lower, yend = upper), size=0.5) +
  geom_hline(yintercept =0, size=1, linetype = "dashed") +
  #xlim(-3,5.5) + 
  ylab("Estimated Change in Log(Poltical Killings)") +
  xlab("UP Vote Share Range") +
  theme_few() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, vjust = 1.1, hjust=1)) +
  scale_color_manual(values=c("black","black","black","black","black","black","black","black","black","black",
                              "black", "black", "black", "black", "black", "black","black"))

mass.cfplot6

#figure 4
ggsave("plot-windows-1.jpg", plot = mass.cfplot6,
       scale = 3, width = 550, height = 504, units = "px",
       dpi = 300, limitsize = TRUE)

#non-dissident victims
c_model20d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= -5 & margin >= -11,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model21d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= -4 & margin >= -10,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model22d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= -3 & margin >= -9,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model23d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= -2 & margin >= -8,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model24d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= -1 & margin >= -7,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model25d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= 0 & margin >= -6,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model26d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= 1 & margin >= -5,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model27d <- lm_robust(log(nonmilitant+1) ~  ifelse(margin <= 2 & margin >= -4,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model28d <- lm_robust(log(nonmilitant+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model29d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= 4 & margin >= -2,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model30d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= 5 & margin >= -1,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model31d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= 6 & margin >= 0,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model32d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= 7 & margin >= 1,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model33d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= 8 & margin >= 2,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model34d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= 9 & margin >= 3,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model35d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= 10 & margin >= 4,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")
c_model36d <- lm_robust(log(nonmilitant+1) ~ ifelse(margin <= 11 & margin >= 5,1,0) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")

mass.plot7 <- matrix(NA,17,3)

mass.plot7[,1] <- c(summary(c_model20d)$coefficients[2,1] - 1.96*summary(c_model20d)$coefficients[2,2],
                    summary(c_model21d)$coefficients[2,1] - 1.96*summary(c_model21d)$coefficients[2,2], 
                    summary(c_model22d)$coefficients[2,1] - 1.96*summary(c_model22d)$coefficients[2,2], 
                    summary(c_model23d)$coefficients[2,1] - 1.96*summary(c_model23d)$coefficients[2,2],
                    summary(c_model24d)$coefficients[2,1] - 1.96*summary(c_model24d)$coefficients[2,2],
                    summary(c_model25d)$coefficients[2,1] - 1.96*summary(c_model25d)$coefficients[2,2],
                    summary(c_model26d)$coefficients[2,1] - 1.96*summary(c_model26d)$coefficients[2,2],
                    summary(c_model27d)$coefficients[2,1] - 1.96*summary(c_model27d)$coefficients[2,2],
                    summary(c_model28d)$coefficients[2,1] - 1.96*summary(c_model28d)$coefficients[2,2],
                    summary(c_model29d)$coefficients[2,1] - 1.96*summary(c_model29d)$coefficients[2,2],
                    summary(c_model30d)$coefficients[2,1] - 1.96*summary(c_model30d)$coefficients[2,2],
                    summary(c_model31d)$coefficients[2,1] - 1.96*summary(c_model31d)$coefficients[2,2],
                    summary(c_model32d)$coefficients[2,1] - 1.96*summary(c_model32d)$coefficients[2,2],
                    summary(c_model33d)$coefficients[2,1] - 1.96*summary(c_model33d)$coefficients[2,2],
                    summary(c_model34d)$coefficients[2,1] - 1.96*summary(c_model34d)$coefficients[2,2],
                    summary(c_model35d)$coefficients[2,1] - 1.96*summary(c_model35d)$coefficients[2,2],
                    summary(c_model36d)$coefficients[2,1] - 1.96*summary(c_model36d)$coefficients[2,2])
mass.plot7[,2] <- c(summary(c_model20d)$coefficients[2,1],
                    summary(c_model21d)$coefficients[2,1], 
                    summary(c_model22d)$coefficients[2,1], 
                    summary(c_model23d)$coefficients[2,1],
                    summary(c_model24d)$coefficients[2,1],
                    summary(c_model25d)$coefficients[2,1],
                    summary(c_model26d)$coefficients[2,1],
                    summary(c_model27d)$coefficients[2,1],
                    summary(c_model28d)$coefficients[2,1],
                    summary(c_model29d)$coefficients[2,1],
                    summary(c_model30d)$coefficients[2,1],
                    summary(c_model31d)$coefficients[2,1],
                    summary(c_model32d)$coefficients[2,1],
                    summary(c_model33d)$coefficients[2,1],
                    summary(c_model34d)$coefficients[2,1],
                    summary(c_model35d)$coefficients[2,1],
                    summary(c_model36d)$coefficients[2,1])
mass.plot7[,3] <- c(summary(c_model20d)$coefficients[2,1] + 1.96*summary(c_model20d)$coefficients[2,2],
                    summary(c_model21d)$coefficients[2,1] + 1.96*summary(c_model21d)$coefficients[2,2], 
                    summary(c_model22d)$coefficients[2,1] + 1.96*summary(c_model22d)$coefficients[2,2], 
                    summary(c_model23d)$coefficients[2,1] + 1.96*summary(c_model23d)$coefficients[2,2],
                    summary(c_model24d)$coefficients[2,1] + 1.96*summary(c_model24d)$coefficients[2,2],
                    summary(c_model25d)$coefficients[2,1] + 1.96*summary(c_model25d)$coefficients[2,2],
                    summary(c_model26d)$coefficients[2,1] + 1.96*summary(c_model26d)$coefficients[2,2],
                    summary(c_model27d)$coefficients[2,1] + 1.96*summary(c_model27d)$coefficients[2,2],
                    summary(c_model28d)$coefficients[2,1] + 1.96*summary(c_model28d)$coefficients[2,2],
                    summary(c_model29d)$coefficients[2,1] + 1.96*summary(c_model29d)$coefficients[2,2],
                    summary(c_model30d)$coefficients[2,1] + 1.96*summary(c_model30d)$coefficients[2,2],
                    summary(c_model31d)$coefficients[2,1] + 1.96*summary(c_model31d)$coefficients[2,2],
                    summary(c_model32d)$coefficients[2,1] + 1.96*summary(c_model32d)$coefficients[2,2],
                    summary(c_model33d)$coefficients[2,1] + 1.96*summary(c_model33d)$coefficients[2,2],
                    summary(c_model34d)$coefficients[2,1] + 1.96*summary(c_model34d)$coefficients[2,2],
                    summary(c_model35d)$coefficients[2,1] + 1.96*summary(c_model35d)$coefficients[2,2],
                    summary(c_model36d)$coefficients[2,1] + 1.96*summary(c_model36d)$coefficients[2,2])

mass.plot7 <- as.data.frame(mass.plot7)
mass.plot7$names <- c("39-45","40-46","41-47","42-48","43-49","44-50","45-51","46-52","47-53","48-54",
                      "49-55","50-56","51-57","52-58","53-59","54-60","55-61")
mass.plot7$names <- as.factor(mass.plot7$names)
colnames(mass.plot7) <- c("lower","ptest","upper","names")

mass.cfplot7 <- ggplot(data = mass.plot7, aes(x = names, y = ptest, color = names)) +
  geom_point(size=2.5) +
  geom_segment(aes(x = names, xend = names, y = lower, yend = upper), size=0.5) +
  geom_hline(yintercept =0, size=1, linetype = "dashed") +
  #xlim(-3,5.5) + 
  ylab("Estimated Change in Log(Non-Dissident Killings)") +
  xlab("UP Vote Share Range") +
  theme_few() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, vjust = 1.1, hjust=1)) +
  scale_color_manual(values=c("black","black","black","black","black","black","black","black","black","black",
                              "black", "black", "black", "black", "black", "black","black"))

mass.cfplot7

#figure 4
ggsave("plot-windows-2.jpg", plot = mass.cfplot7,
       scale = 3, width = 550, height = 504, units = "px",
       dpi = 300, limitsize = TRUE)

#3.6 placebo: violent dissidents  ####

c_model7a <- lm_robust(log(militant_violent + 1) ~ close_election_3pt, data=chile.commune,se_type="HC2")

c_model7b <- lm_robust(log(militant_violent + 1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                         as.factor(province) + armybase, data=chile.commune, se_type="HC2")

#compare dissident coefficients to non-dissident outcome coefficients
ztest <- (c_model1d$coefficients[2] - c_model7b$coefficients[2])/
          sqrt(c_model1d$std.error[2]^2 + c_model7b$std.error[2]^2)

#two-sided p-value
2*pnorm(ztest, mean = 0, sd = 1, lower.tail = FALSE)


#3.7 robustness: sensitivity analysis ####
#using cinelli and hazlett (2020)

sens1 <- sensemakr(lm(log(totalvictims+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                        as.factor(province) + armybase, data=chile.commune),
                   treatment = "close_election_3pt", benchmark_covariates = "log(pop_1970_est)")

sens1$sensitivity_stats

sens1ovb <- ovb_bounds(lm(log(totalvictims+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                            as.factor(province) + armybase, data=chile.commune),
                       treatment = "close_election_3pt", benchmark_covariates = "log(pop_1970_est)")

sens1ovb

sens1ovb_plot <-  ovb_contour_plot(lm(log(totalvictims+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                                        as.factor(province) + armybase, data=chile.commune),
                                   treatment = "close_election_3pt", benchmark_covariates = "log(pop_1970_est)")

sens2 <- sensemakr(lm(log(nonmilitant+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                        as.factor(province) + armybase, data=chile.commune),
                   treatment = "close_election_3pt", benchmark_covariates = "log(pop_1970_est)")

sens2$sensitivity_stats

sens2ovb <- ovb_bounds(lm(log(nonmilitant+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                            as.factor(province) + armybase, data=chile.commune),
                       treatment = "close_election_3pt", benchmark_covariates = "log(pop_1970_est)")

sens2ovb

sens2ovb_plot <-  ovb_contour_plot(lm(log(nonmilitant+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                                        as.factor(province) + armybase, data=chile.commune),
                                   treatment = "close_election_3pt", benchmark_covariates = "log(pop_1970_est)")


#3.8 mechanism: participation in repression ####

chile.commune$rettig_collab <- ifelse(chile.commune$comuna %in% c("salamanca","santa barbara","quilaco","quilleco",
                                                                  "mulchen","liquire","entre lagos","paine"),1,0) 

mean(chile.commune$up_1973_pct[chile.commune$rettig_collab == 1], na.rm=T)
mean(chile.commune$up_1973_pct[chile.commune$rettig_collab == 0], na.rm=T)

t.test(x = chile.commune$up_1973_pct[chile.commune$rettig_collab == 1], 
       y = chile.commune$up_1973_pct[chile.commune$rettig_collab == 0], alternative = "greater", mu = 0, 
       paired = FALSE, var.equal = FALSE, conf.level = 0.95)

t.test(x = chile.commune$nonmilitant[chile.commune$rettig_collab == 1], 
       y = chile.commune$nonmilitant[chile.commune$rettig_collab == 0], alternative = "greater", mu = 0, 
       paired = FALSE, var.equal = FALSE, conf.level = 0.95)

t.test(x = chile.commune$totalvictims[chile.commune$rettig_collab == 1], 
       y = chile.commune$totalvictims[chile.commune$rettig_collab == 0], alternative = "greater", mu = 0, 
       paired = FALSE, var.equal = FALSE, conf.level = 0.95)

ks.test(x = chile.commune$up_1973_pct[chile.commune$rettig_collab == 1], 
        y = chile.commune$up_1973_pct[chile.commune$rettig_collab == 0], alternative = "less",
        exact = NULL)

chile.commune$rettig_collab <- as.factor(chile.commune$rettig_collab)
levels(chile.commune$rettig_collab) <- c("No Participation","Participation")

rettig_plot <- ggplot(chile.commune, aes(up_1973_pct, stat(density), linetype = rettig_collab)) +
  geom_freqpoly(binwidth = 10,lwd=1.25) +
  ylab("Density") +
  xlab("UP 1973 Vote Share") +
  scale_linetype_discrete(name = "Civilian Participation") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"),
        legend.position="bottom")

rettig_plot

#figure 5
ggsave("participation.jpg", plot = rettig_plot,
       scale = 3, width = 600, height = 600, units = "px",
       dpi = 300, limitsize = TRUE)

#mediation analysis

chile.commune$rettig_collab <- ifelse(chile.commune$comuna %in% c("salamanca","santa barbara","quilaco","quilleco",
                                                                  "mulchen","liquire","entre lagos","paine"),1,0) 

c_model_mediate1 <- lm(rettig_collab ~ abs(margin) , data=chile.commune)
c_model_mediate2 <- lm(log(totalvictims+1) ~ abs(margin) + rettig_collab, data=chile.commune)

mediation.model1 <- mediate(c_model_mediate1, c_model_mediate2,treat = "abs(margin)",
                           mediator = "rettig_collab",sims=10000)
summary(mediation.model1)

mediation.sens <- medsens(mediation.model1, rho.by = 0.1, effect.type = "indirect", sims = 100)
plot(mediation.sens,ylim=c(-0.02,0.02))

#3.9 robustness: balcells index ####

chile.commune$balcells <- 1 - ((chile.commune$up_1973_pct - (100 - chile.commune$up_1973_pct))/100)^2

#total victims outcome
c_model41a <- lm_robust(log(totalvictims+1) ~ ifelse(is.infinite(log(1-balcells)),NA,-log(1-balcells)), data=chile.commune,se_type="HC2")

c_model41b <- lm_robust(log(totalvictims+1) ~ ifelse(is.infinite(log(1-balcells)),NA,-log(1-balcells)) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")

#non-militant outcome
c_model41c <- lm_robust(log(nonmilitant+1) ~ ifelse(is.infinite(log(1-balcells)),NA,-log(1-balcells)), data=chile.commune,se_type="HC2")

c_model41d <- lm_robust(log(nonmilitant+1) ~ ifelse(is.infinite(log(1-balcells)),NA,-log(1-balcells)) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")

#3.10 robustness: continuous vote margin ####

#total victims outcome

c_model42a <- lm_robust(log(totalvictims+1) ~ abs(margin), data=chile.commune,se_type="HC2")

c_model42b <- lm_robust(log(totalvictims+1) ~ abs(margin) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")

#non-militant outcome
c_model42c <- lm_robust(log(nonmilitant+1) ~ abs(margin), data=chile.commune,se_type="HC2")

c_model42d <- lm_robust(log(nonmilitant+1) ~ abs(margin) + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")

#3.11 robustness: varying definition of militancy ####

#strict (only leaders (political and labor) and violent nonstate groups removed)

c_model46a <- lm_robust(log(totalvictims - militant_strict + 1) ~ close_election_3pt, data=chile.commune,se_type="HC2")

c_model46b <- lm_robust(log(totalvictims - militant_strict + 1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")

#loose (remove students too)

c_model46c <- lm_robust(log(totalvictims - militant_loose + 1) ~ close_election_3pt, data=chile.commune,se_type="HC2")

c_model46d <- lm_robust(log(totalvictims - militant_loose + 1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")


#3.12 alternative: caravan of death ####

#caravan of death

#militant outcome
c_model38a <- lm_robust(log(totalvictims - nonmilitant +1) ~ caravana_muerte, data=chile.commune,se_type="HC2")

c_model38b <- lm_robust(log(totalvictims - nonmilitant +1) ~ caravana_muerte + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")

#non-militant outcome
c_model38c <- lm_robust(log(nonmilitant+1) ~ caravana_muerte, data=chile.commune,se_type="HC2")

c_model38d <- lm_robust(log(nonmilitant+1) ~ caravana_muerte + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")

#count
#militant outcome
c_model40a <- glm.nb(totalvictims - nonmilitant ~ caravana_muerte, data=chile.commune)

c_model40b <- glm.nb(totalvictims - nonmilitant ~ caravana_muerte + allende_1970_pct + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                       as.factor(province) + armybase, data=chile.commune)

#nonmilitant outcome
c_model40c <- glm.nb(nonmilitant ~ caravana_muerte, data=chile.commune)

c_model40d <- glm.nb(nonmilitant ~ caravana_muerte + allende_1970_pct + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                       as.factor(province) + armybase, data=chile.commune)


#3.13 alternative: without turned-in dissidents ####

chile.commune$total_victims_alt <- chile.commune$totalvictims - chile.commune$militant_turnedin

c_model48a <- lm_robust(log(total_victims_alt+1) ~ close_election_3pt, data=chile.commune,se_type="HC2")

c_model48b <- lm_robust(log(total_victims_alt+1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=chile.commune, se_type="HC2")

#3.14 alternative: without large comunas ####

#without capital

c_model49a <- lm_robust(log(totalvictims + 1) ~ close_election_3pt, data=subset(chile.commune,chile.commune$capital==0),se_type="HC2")

c_model49b <- lm_robust(log(totalvictims + 1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=subset(chile.commune,chile.commune$capital==0), se_type="HC2")

c_model49c <- lm_robust(log(nonmilitant + 1) ~ close_election_3pt, data=subset(chile.commune,chile.commune$capital==0),se_type="HC2")

c_model49d <- lm_robust(log(nonmilitant + 1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=subset(chile.commune,chile.commune$capital==0), se_type="HC2")

#without urban areas 

c_model50a <- lm_robust(log(totalvictims + 1) ~ close_election_3pt, data=subset(chile.commune, chile.commune$pop_1970_est < 100000),se_type="HC2")

c_model50b <- lm_robust(log(totalvictims + 1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=subset(chile.commune,chile.commune$pop_1970_est < 100000), se_type="HC2")

c_model50c <- lm_robust(log(nonmilitant + 1) ~ close_election_3pt, data=subset(chile.commune,chile.commune$pop_1970_est < 100000),se_type="HC2")

c_model50d <- lm_robust(log(nonmilitant + 1) ~ close_election_3pt + log(allende_1970_votes) + polarized + infant_mortality_pct + log(pop_1970_est) + votershare +
                          as.factor(province) + armybase, data=subset(chile.commune,chile.commune$pop_1970_est < 100000), se_type="HC2")

#4 tables####

#4.1 descriptive tables####

desc.stats <- function(var){
  stats <- vector(mode="numeric")
  stats[1] <- round(min(var, na.rm=T),2)
  stats[2] <- round(max(var, na.rm=T),2)
  stats[3] <- round(mean(var, na.rm=T),2)
  stats[4] <- round(median(var, na.rm=T),2)
  stats[5] <- round(sd(var, na.rm=T),2)
  stats[6] <- sum(ifelse(is.na(var),1,0))
  return(stats)}

desc.covs_commune <- as.data.frame(
  cbind(c("Total Victims","Non-Militant Victims","Close Election, 1973", "Log(1970 Population)",
          "Infant Mortality, 1972","Allende Support","Polarized","Voter Share",
          "Caravan of Death","Army Base","Log(Marriages, 1972)","Catholic Share, 1968"),
        rbind(desc.stats(chile.commune$totalvictims),
              desc.stats(chile.commune$nonmilitant),
              desc.stats(chile.commune$close_election_3pt),
              desc.stats(log(chile.commune$pop_1970_est)),
              desc.stats(chile.commune$infant_mortality_pct),
              desc.stats(chile.commune$allende_1970_pct),
              desc.stats(chile.commune$polarized),
              desc.stats(chile.commune$votershare),
              desc.stats(chile.commune$caravana_muerte),
              desc.stats(chile.commune$armybase),
              desc.stats(log(chile.commune$marriages+1)),
              desc.stats(chile.commune$catholicity_1968))))

colnames(desc.covs_commune) <- c("Variable","Minimum","Maximum","Mean","Median","Std. Dev","Missing")

print(xtable(desc.covs_commune, caption="Descriptive Statistics, Municipality Level"), include.rownames=FALSE) #UPDATED
print(xtable(bal, caption="Comparing Covariates, Competitive and Uncompetitive"), include.rownames=FALSE) #UPDATED


#4.2 balance table####

balance <- function(samp.bal,treat){
  treatnum <- nrow(filter(samp.bal, treat == 1))
  controlnum <- nrow(filter(samp.bal, treat == 0))
  
  names <- colnames(samp.bal)
  n <- ncol(samp.bal)
  diff_mean <- array(NA,c(1,n))
  se <- array(NA,c(1,n))
  t <- array(NA,c(1,n))
  p <- array(NA,c(1,n))
  
  for (i in 1:n){
    diff_mean[i] <- mean(filter(samp.bal, treat == 1)[,i],na.rm=T) - 
      mean(filter(samp.bal, treat == 0)[,i],na.rm=T)
    se[i] <- sqrt(var(filter(samp.bal, treat == 1)[,i],na.rm=T)/treatnum +
                    var(filter(samp.bal, treat == 0)[,i],na.rm=T)/controlnum)
    t[i] <- diff_mean[i]/se[i]
    p[i] <- 2*pt(-abs(t[i]), df=pmin(treatnum,controlnum)-1)
  }
  t.bal <- t(rbind(names,diff_mean,t,se,p))
  colnames(t.bal) <- c("variable","diff_mean","t","se","p")
  t.bal[,-1] <- round(as.numeric(t.bal[,-1]),4)
  t.bal
}

dat.balance <- as.data.frame(cbind(chile.commune$close_election_3pt,chile.commune$votershare, chile.commune$allende_1970_pct,
                                   chile.commune$armybase, chile.commune$infant_mortality_pct,
                                   log(chile.commune$pop_1970_est), log(chile.commune$marriages+1)))

colnames(dat.balance) <- c("close_election_3pt","votershare","allende_1970_pct","armybase","infant_mortality_pct",
                           "log(pop_1970_est)","log(marriages)")

bal <- as.data.frame(balance(samp.bal = dat.balance, treat = dat.balance$close_election_3pt))

#4.3 regression tables####

texreg(list(c_model1a,c_model1b,c_model1c,c_model1d), include.ci = FALSE, omit.coef = "as.factor") #main results - table 1

#online appendix tables
texreg(list(c_model3a,c_model3b,c_model3c,c_model3d), include.ci = FALSE, omit.coef = "as.factor") #count models 
texreg(list(c_model4a,c_model4b,c_model4c,c_model4d), include.ci = FALSE, omit.coef = "as.factor") #ZINB models 
texreg(list(c_model5a,c_model5b,c_model5c,c_model5d), include.ci = FALSE, omit.coef = "as.factor") #killing rate models 
texreg(list(c_model13a,c_model13b, c_model13c, c_model13d), include.ci = FALSE, omit.coef = "as.factor") #september 13 and later 
texreg(list(c_model35c,c_model34c,c_model33c,c_model32c,c_model31c,c_model30c,c_model29c,c_model28c), include.ci = FALSE, omit.coef = "as.factor") #shifting thresholds 1 
texreg(list(c_model27c,c_model26c,c_model25c,c_model24c,c_model23c,c_model22c,c_model21c,c_model20c), include.ci = FALSE, omit.coef = "as.factor") #shifting thresholds 2 
texreg(list(c_model35d,c_model34d,c_model33d,c_model32d,c_model31d,c_model30d,c_model29d,c_model28d), include.ci = FALSE, omit.coef = "as.factor") #shifting thresholds 3 
texreg(list(c_model27d,c_model26d,c_model25d,c_model24d,c_model23d,c_model22d,c_model21d,c_model20d), include.ci = FALSE, omit.coef = "as.factor") #shifting thresholds 4 
texreg(list(c_model35,c_model34,c_model33,c_model32,c_model31,c_model30,c_model29,c_model28), include.ci = FALSE, omit.coef = "as.factor") #widening thresholds 1 
texreg(list(c_model27,c_model26,c_model25,c_model24,c_model23,c_model22,c_model21,c_model20), include.ci = FALSE, omit.coef = "as.factor") #widening thresholds 2 
texreg(list(c_model35a,c_model34a,c_model33a,c_model32a,c_model31a,c_model30a,c_model29a,c_model28a), include.ci = FALSE, omit.coef = "as.factor") #widening thresholds 3 
texreg(list(c_model27a,c_model26a,c_model25a,c_model24a,c_model23a,c_model22a,c_model21a,c_model20a), include.ci = FALSE, omit.coef = "as.factor") #widening thresholds 4 
texreg(list(c_model38a,c_model38b,c_model38c,c_model38d), include.ci = FALSE, omit.coef = "as.factor") #caravan of death 
texreg(list(c_model40a,c_model40b,c_model40c,c_model40d), include.ci = FALSE, omit.coef = "as.factor") #caravan of death, count 
texreg(list(c_model41a,c_model41b,c_model41c,c_model41d), include.ci = FALSE, omit.coef = "as.factor") #balcells measurement 
texreg(list(c_model42a,c_model42b,c_model42c,c_model42d), include.ci = FALSE, omit.coef = "as.factor") #margin of victory 
texreg(list(c_model7a,c_model7b), include.ci = FALSE, omit.coef = "as.factor") #placebo outcome - violent dissidents 
texreg(list(c_model46a,c_model46b,c_model46c,c_model46d), include.ci = FALSE, omit.coef = "as.factor") #changing dissident classification 
texreg(list(c_model48a,c_model48b), include.ci = FALSE, omit.coef = "as.factor") #exclude turned in dissidents 
texreg(list(c_model49a,c_model49b,c_model49c,c_model49d), include.ci = FALSE, omit.coef = "as.factor") #excluding santiago 
texreg(list(c_model50a,c_model50b,c_model50c,c_model50d), include.ci = FALSE, omit.coef = "as.factor") #excluding urban areas over 100,000 population 
texreg(list(c_model13e,c_model13f,c_model13g,c_model13h,c_model13i,c_model13j,c_model13k,c_model13l), include.ci = FALSE, omit.coef = "as.factor") #effects by month 
